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We study the quantum error correction tiiresiiold of Kitaev's toric code over the group Zd subject 
to a generalized bit-flip noise. This problem requires novel decoding techniques, and for this purpose 
we generalize the renormalization group method we previously introduced in [TJ[2] for Z2 topological 
codes. 
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Kitaev's topological code (KTC) [5] on qubits is the 
archetypical topological code and has been extensively 
studied. As explained in Kitaev's original paper [3|, this 
construction applies to any group. Much less is known 
about these generalizations, and in this paper we inves- 
tigate the quantum error correction (QEC) thresholds of 
the KTCs built with the groups Z^, where d > 2. We 
label these as Z^j-KTC, so the original code on qubits 
corresponds to Z2-KTC. 

As explained in [4 , Z2-KTC can be decoded by a bi- 
nary perfect matching algorithm [5^, since every particle 
is its own anti-particle in this model. Because this is 
not the case for d > 2, other techniques are required 
and for this purpose we generalize the renormalization 
group (RG) soft decoder that we introduced in [TJ [2]. 
Our numerical simulations show that the threshold in- 
creases monotonically with d and appears to follow the 
general trend of the qudit hashing bound. 

This paper is organized as follows. First, we introduce 
a generalized Pauli group (see [SI [7] for more details), 
stabilizer codes, and Z^-Kitaev's toric code. Next, we 
briefly review the decoding problem of these systems and 
show how the RG decoder applies in this case. Finally, we 
present the numerical results and close with a discussion. 



I. GENERALIZATION OF KITAEV'S TORIC 
CODE 

In this section, we review the definition of Z^-KTC 
and show that many features of KTC on qubits extend to 
them. Since we will be working with qudits, we introduce 
a generalized Pauli group. The Hilbert space of a qudit, 
Hd, is spanned by the states {|0), |1), . . . , |d — 1)}. We 
define the operators X and Z such that 

^|5> = I.9®1), (1) 
Z\9)=uj^\9), 

where < g < d, "0" denotes addition modulo d, and 
to = qI-'^t'/<i_ Tiig generalized Pauli group is generated by 
X, Z, and a phase, i.e., Vd = {^,X, Z) if d is odd and 
-p^ = (cj^/^ A:,Z) if d is even {XZ has order 2d in this 
case). From the definitions of Eq. (fTl), we deduce the 
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FIG. 1: Zd-KTC stabilizer generators. To each vertex v, we 
associate an operator A„(left) and to each plaquette p, we 
associate an operator Bp (right). 

following properties 

^1ff) = l5®a), 

Z'^\g)=uj'^s\g), (2) 
ZX\g)^iuXZ\g), 
Z°'X^\g) ^uj°'^X^Z''\g). 

Lastly, we define the n-qudit Pauli group V^^ = Vf" as 
the n-fold tensor product of Vd- 

The stabilizer group S is an ablian subgroup of 1^2- 
The code is defined as the simultaneous +1 eigenspace 
of all stabilizers. Note that even though the general- 
ized Pauli operators are unitary, they are not hermitian 
in general so do not correspond to physical observables. 
However, the operator |(s -I- s^) is hermitian and can be 
measured. Since s has eigenvalues w", ^(s-|-s^) has eigen- 
values ^{u)"" +uj~°') = cos(27ra/(i) which are in one-to-one 
correspondence with the eigenvalues of s. 

With these definitions in place, we present a general- 
ization of KTC on qudits, which we call Z^j-KTC, using 
Kitaev's original construction [3] on the cyclic groups Zd 
with d > 2. The system is a square lattice of linear size 
L with periodic boundary conditions. Each edge is occu- 
pied by a qudit, so there are in total n = 2L^ qudits. We 
define vertex operators Ay and plaquette operators Bp 
as shown in Fig. [T] There is one such operator for each 
vertex and each plaquette. We verify that they commute 
using the last line of Eq. ([2]). These operators generate 
the stabilizer group S = {A^, Bp) and the code is spanned 
by the simultaneous +1 eigenstates of the stabilizer gen- 
erators. 

Figure [2] illustrates how applying some power of X on 
a codestate creates defects on the lattice. Indeed, 
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FIG. 2: Plaquette defects created by the application of some 
power of X. The values a {—a) in the plaquettes are such 
that the eigenvalue of the corresponding Bp is oj" By 
choosing appropriately the powers of X, we can build string 
operators with defects only on their endpoints. Non-trivial 
cocyles of correspond to X logical operators. 



applied on some qudit does not commute with the two 
plaquette operators involving that qudit. The eigenval- 
ues of the plaquettes to the north or east of the error will 
change from 1 to and those of the plaquettes to the 
south or west will change from 1 to One can show 

that the defects thus created are topological charges; we 
associate the charge a to a plaquette defect corresponding 
to an eigenvalue of that plaquette. With this choice 
of labeling, the charge group restricted to plaquettes is 
Zd with addition. 

From these simple facts, it follows that string operators 
can be built with defects attached only to their endpoints 
(these strings actually live on the dual lattice, just like 
in KTC). This requires a careful choice of the powers 
of X on the qudits along the string such that the total 
charge in each plaquette is except on its endpoints. For 
instance, one can adopt the convention that power a is 
used when heading north or east, and —a when heading 
south or west. Moreover, we can verify that non-trivial 
cocycles (loops on the dual lattice, see Fig. [2]) of any 
power of X obeying this convention commute with the 
stabilizer. These operators are not in the stabilizer as 
all the vertex generators of Fig. [T] are trivial cocyles. It 
follows that such operators, e.g. the one found at the 
bottom of Fig. [2j are logical operators (for any value of 
a). 

A similar analysis holds for defects created by powers of 
Z operators. In this case, the defects live on vertices and 
string operators, on the direct lattice. Also, non-trivial 
cycles of any power of Z are logical operators. From 
the form of the logical operators, we directly deduce that 
there are two qudits encoded in the code space. Again, 
this is analogous to the case of KTC. 



II. Zd-KTC DECODING 



We are now interested in the problem of error cor- 
recting Zrf-KTCs for d > 2. In our study, we consider 
a simple noise model that generalizes the independent 
symmetric bit-flip channel to qudits^: with probability 
1 — Pphys J the qudit remains unaffected and with proba- 
bility pphys, we apply at random (uniformely distributed) 
one of X, X"^, . . . , X'^~^. Suppose an error E e oc- 
curs on a code state. It creates defects on the lattice 
and by measuring the eigenvalues of every ^{Ay + Al) 
and ^{Bp + Bj^) we can learn the position and charge 
of each defect. The role of the decoder is to bring the 
system back in the code space by applying a correcting 
Pauli operator, C € V^- However, care must be taken in 
choosing an appropriate correcting operation. Indeed, if 
the operator CE resulting from the combination of the 
error and the recovery is an element of S, the state is 
unaffected. However, if CE is a non-trivial logical oper- 
ator, then the system is returned to the code space but 
potentially in a different code state, so the information 
is corrupted. 



Any operator E £ creating the measured configu- 
ration of defects is a potential error. However, we classify 
these operators by their logical effect on the code space: 
two operators Ei , E2 with the same configuration of de- 
fects are equivalent iff e\Ei has a trivial effect on the 
code, i.e. Ei ^ E2 iff e\Ei G S. Note that since Ei and 
E2 lead to the same defect configuration, E\Ei creates 
no defect, or equivalently. Ex creates some defects that 
e\ annihilates. 



Given a measured defect configuration, the decoder 
seeks for the best correction among the set of all errors 
which would lead to this defect configuration. One strat- 
egy would be to identify the error from this set that has 
the largest probability V{E), where the probability of an 
error is specified by the physical noise model, in our case 
the symmetric bit-flip channel. This turns out not to be 
optimal however, because some errors have equivalent ef- 
fects on all code states. Thus, the decoder should instead 
seek for the most likely equivalence class of errors. The 
probability of an equivalence class of errors is obtained 
by summing over the probability of each error within a 
class. Given these probabilities, the optimal correction 
consists in applying the adjoint of any representative of 
the class with maximal probability. 



^ This noise model can also be seen as emerging from a qudit 
depolarization channel that maps p— > {l — q)p + q^ when X and 
Z errors are treated independently, and Pphys = d^^). 
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FIG. 3: (a) The lattice is cut into unit cells containing ten 
qudits (edges). The renormalization process takes the defect 
configuration and the noise model on a unit cell as inputs 
and outputs a two-qudit distribution (white disks) which cor- 
responds to a probability on the charge flow through the cor- 
responding boundaries. Green disks represent plaquette op- 
erators. The plaquette corresponding to the green circle is 
replaced by the product of all four plaquettes of the unit cell, 
such that its eigenvalue gives the total charge of the cell. This 
value is only going to be used in the next round of RG (larger 
green disk), (b) Labeling convention for qudits in Eq. (|3| 



III. RG DECODER GENERALIZATION TO 

Zd-KTC 

Unfortunately, the above procedure cannot be realized 
efficiently in general since the number of errors in each 
equivalence class scales exponentially with the system 
size. In [HE], we introduced a renormalization group 
soft decoder (RG decoder) that efficiently approximates 
the exact calculation (see ^Sj for a related scheme). The 
general idea is to cut the lattice into small unit cells (e.g. 
2x2 sub-lattices) and to "distill" from each cell an effec- 
tive two-qubit noise model, c.f. Fig.jsj^a). This is realized 
by keeping track of the flow of charges through the cell 
and summing over the microscopic details leading to this 
flow. This has the effect of shrinking the lattice linear 
size by a constant factor (fc for cells of size k x k). Re- 
cursing on this process, one can shrink the lattice to a 
constant, manageable, size where the exact decoding can 
be performed. With appropriate simple modifications, 
this method can be used for charges over Z^. 

There are two technical difficulties in realizing the 
above heuristic description, which are both caused by 
charge conservation. First, because the unit cells share 
boundaries, the flow of charge through one boundary of 
a cell should be equal and opposite to the flow of charge 
of the corresponding boundary of the neighbouring cell. 
Thus, the variable corresponding to charge flows in each 
cell are highly constrained. This problem is easily cir- 
cumvented by keeping only track of the flow of charge 
through the northern and the western boundary of each 
cell, i.e. by eliminating this redundancy. 

Second, the sum of the charge flow through the bound- 
aries of a cell must be equal to its total charge, revealed 
by the syndrome measurement. This once again sets a 
hard constraint between the variables corresponding to 
the charge flows, which would in principle require a prob- 
ability distribution that correlates all the variables of the 



system. This cannot be realized efficiently, so we must re- 
sort to some approximation. As a first approximation, we 
choose to ignore the cross-cell correlations, and keep only 
marginal probabilities on the flows associated to a given 
cell (we keep a probability distribution that involves the 
northern and western boundary only). To diminish the 
effect of these correlations we are neglecting, we let the 
charge inside a unit cell fluctuate. For each unit cell, we 
measure all but one of the plaquettes it encloses. This re- 
maining plaquette thus determines the total charge of the 
unit cell, and indeed we can substitute the corresponding 
stabilizer generator by a plaquette enclosing the entire 
unit cell (obtained by multiplying all the plaquette op- 
erators contained in the unit cell). This new stabilizer 
generator represents a renormalized charge. 

This procedure is illustrated on Fig. [sj^a) where green 
disks represent plaquettes that are measured and the 
green circle represents the plaquette that is left fluctu- 
ating. This green circle is replaced by the larger, renor- 
malized green disk (on the right) that is used in the next 
RG step. The white disks on this figure each represent a 
probability distribution on charge flow, or equivalently a 
two-qudit probability distribution. Thus, after one round 
of RG, we are left with a smaller lattice and both renor- 
malized charges and renormalized noise models. 

Equation ([s]) lists a set of generators for all X oper- 
ators living on a unit cell (see Fig. [sf^b) for labelling). 
This basis will be used to decompose any X-type error 
contained on the unit cell. These operators are defined 
in accordance to the renormalization process itself as we 
now explain. The operators are used to build a repre- 
sentative error with the appropriate defect configuration. 
Indeed, only the Ti operators of Eq. ^ do not commute 
with all three plaquette operators in the unit cell (green 
disks of Fig. [sj^a)). Label the defect configuration on a 
unit cell as a = (oq, ai, 02), where uq is the charge of the 
north-west plaquette, ai is the charge of the north-east 
one, and 02 is the charge of the south-west one. Then, 
the Pauh operator t{a) = T^°T^^T!^^ creates the defect 
configuration a. Moreover, given a defect configuration 
a, every potential error has to contain this product in its 
decomposition on basis Eq. ([S]) since only the opera- 
tors do not commute with plaquettes. The Li operators 
characterize the flow of charge through the northen and 
western boundaries, so the two-qudit ouput distribution 
of a RG round is precisely the probability distribution 
over these two operators. The Si operators are stabilizer 
operators (or parts of stabilizer generators supported on 
the unit cell). They only deform strings without chang- 
ing their defect configuration or their associated charge 
flow. Lastly, the Ei operators correspond to charge flow- 
ing through the southern and eastern boundaries into 
the plaquette operator that is left out. Thus, they are 
responsible for the charge fluctuation inside the unit cell 
and they are summed over. 
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With these definitions, we can formahy describe a RG 
round that starts with a defect configuration a, and com- 
putes the marginal probabihty of each I e {Lq,Li) con- 
ditioned on the measured defect configuration, 



(4) 



where i = T^oT^^T"^ is given by the defect configura- 
tion and Titles) is the probabihty assigned to the error 
E = ties by the noise modeL The complexity of decoding 
a unit cell is given by the number of operators that are 
considered in Eq. \{Lo,Li)\-\{Eo,Ei)\-\{So,Si,S2)\. 
Since all Li, Ei and Si have order d, the complexity is 
the constant . For different unit cell sizes, the com- 
plexity is still a power of d, but with a different exponent 
which depends on the number of qudits in the cell and 
the number of measured stabilizer generators. Moreover, 
the number of unit cells to decode in a given round of 
RG is given by (L/fc)^ where k and L are the linear sizes 
of the unit cell and the global lattice, respectively. Thus, 
the complexity of a step of RG goes as d'^{L/k)'^ for some 
constants c and k that depend on the choice of unit cell. 
Of course, the RG calculations on different cells can be 
executed in parallel. 

The procedure we have described above to evade the 
correlations caused by local charge conservation is only a 
heuristic, and can be improved using belief propagation 
(BP). Roughly, the role of BP is to ensure consistency be- 
tween the marginal probability of qubits located at the 
boundary of two or more unit cells, e.g. qudits 0, 1, 8 and 
9 (see Fig.[3](b) for labeling). First, given a defect config- 
uration inside a unit cell, one can compute the marginal 
error probability Vq{tles\q) for each qudit q, obtained by 
taking a marginal of Titles). These are called messages 
and denoted m°^^{p), where q labels a qudit and p is a 
one-qudit Pauli operator. These outgoing messages are 
then exchanged between neighbouring cells, and become 
incoming messages, e.g. a cell c sends to its northern 
neighbour c' the message mg"* that becomes mjf in c', 
and receives from c' the message mg"* that becomes TOq" 
in c. Subsequent rounds of messages can be calculated 
using the received messages, following the prescription 
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FIG. 4: Threshold estimation for Z3-KTC. The x-axis repre- 
sents physical error rate and the y-axis, decoding error rate. 
The blue dots, red squares and yellow diamonds correspond 
to L = 32, L — 64 and L = 128 respectively. The fitting 
curve used is pdec = ipphys ~~ pth)L^^" ■ In this case, we find 
Pth = 0.13(0). 
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FIG. 5; The blue diamonds are the values extracted by fit- 
ting the threshold values for 2 < n < 6 (see Fig. |4] for 
example). The red squares are obtained via the general- 
ized hashing bound (see text) rescaled by a common factor 
a = Pth(2)/C2 ~ 0.81. The error bars are (pessimistically) 
obtained e.g. by replacing each line in Fig. |4] by a stripe of 
width equal to the statistical error bars, and determining the 
values of Pphys above and below the crossing point where the 
strips cease to overlap. We do not report the fitting parame- 
ter V because they are too sensitive to statistical fluctuations 
and therefore unreliable in our study. 



Here, q,q' G {0,1,8,9}, tles\q is the restriction to qudit 
q of the Pauli operator ties and Tq is the marginal on 
qudit q of the noise model as above. BP can be iterated 
a few times (e.g. three rounds) before executing a RG 
step. This has the effect of replacing Eq. Q by 



(5) 
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Titles)Y[m]^itles\q). (6) 
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IV. NUMERICAL RESULTS 

In this section, we present our numerical estimates of 
the thresholds of Z^-KTCs for 2 < d < 6 subject to the 
generalized bit-flip noise model introduced in the pre- 
vious section. The threshold is defined as the value of 
the physical noise rate Pphys below which the decoding 
error probability p^ec can be made arbitrarily small by 
increasing the lattice size L. 

The simulations were performed as follows. For various 
values of d, L and Pphys, specifiying a "Zd-^^G of linear 
size L subject to a noise of parameter Pphys, we performed 
a Monte Carlo simulation to estimate the decoding error 
probability pdcc- We used sample sizes of the order of 
10*. For a fixed value of d, we plotted estimates of pdcc 
vs Pphys for different values of L. We then used the fitting 
model pdcc = (Pphys— Pth)^^^'^ (see [1112] for more details) 
to estimate the value of the threshold. As an example, 
we plotted the results and the fits for Z3-KTC on Fig. 4] 

Repeating this for 3 < d < 6 (2 was studied in [D [2 ), 
Fig. [5] shows pth as a function of d. Heuristically, we did 
expect that the value of pth increases with d. Indeed, 
if we imagine simulating a qudit using log2 d qubits, a 
fixed noise rate for increasing values of d translates into a 
decreased noise rate per qubit. Moreover, it was reported 
in [TU] that the performance of BP for Z^-KTC, which 
is very poor in the qubit case, is greatly increased as d 
grows. 

It is intringuing to note that for Z2-KTC subject to bit- 
flip or depolarizing noise, pth is numerically very close 
to the hashing bound [H [9l |TT]. The hashing bound, 
obtained by a simple packing argument 12J, states that 
for non-degenerate CSS codes, 



0< l-2i?2(p), 



(7) 



where H2 is the binary entropy: i?2(p) = (1 — p) log2(l — 
p) + p log2 p- From Eq. ([t]) , one can calculate the satu- 
rating point C2 ~ 0.110 which is indeed quite close to 
the optimal threshold of the Z2-KTC subject to inde- 
pendent bit-flip and phase-flip errors, pt/i(2) w 0.109(4) 
[3112]. This near coincidence is intriguing given that topo- 
logical codes are highly degenerate, so there is no rea- 
son they should obey the hashing bound. Of course, 
the decoder we are using here is sub-optimal, so the 
threshold we find pth{2) ~ 0.89(6) is a smaller fraction 
a = pth(2)/C2 ~ 0.81(4) of the hashing bound. 
For qudits, the hashing bound is 



In this case, we find C3 « 0.159,6*4 ~ 0.189 and so on. 
Figure[5]shows the threshold pth(d) obtained with the RG 
decoder as well as a rescaled hashing bound aCd where 
a is determined by the Z2 fit. The agreement is both 
unexplained and surprisingly good. Note also that even 
though our decoder is sub-optimal, pth{d + 1) > Cd for 
all d we have studied, which strongly support the claim 
that the threshold increases with d. 



Q<l-2Hd{p), 
with Hd{p) = (1 - p) log(l -p)+p log 



(8) 



d-l 



V. CONCLUSION 

In this paper, we presented a generalization of the 
renormalization group decoder of [U [2] to Kitaev topo- 
logical codes built with the groups Z^. Our numerical 
results show that the threshold value increases as a func- 
tion of the local dimension d. Moreover, its behaviour is 
in very good agreement with a scaling predicted by the 
hashing bound. This trend could be confirmed by more 
accurate numerical estimates using a mapping to a statis- 
tical mechanics model, which does not require solving the 
decoding problem [H [TT] . A theoretical understanding of 
this behavior is also desirable. Lastly, estimating the 
threshold in the presence of measurement error and de- 
tailed syndrome measurement circuits on qudits remains 
an interesting open question. 
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